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Abstract: Both hard and soft QCD dynamics are important in charmonium 
production, as presented here through a next-to-leading order QCD matrix 
element calculation combined with the colour evaporation model. Observed 
xf and p± distributions of J/ip in hadroproduction at fixed target and pp 
collider energies are reproduced. Quite similar results can also be obtained 
in a more phenomenologically useful Monte Carlo event generator where the 
perturbative production of cc pairs is instead obtained through leading order 
matrix elements and the parton shower approximation of the higher order 
processes. The soft dynamics may alternatively be described by the soft colour 
interaction model, originally introduced in connection with rapidity gaps. We 
also discuss the relative rates of different charmonium states and introduce an 
improved model for mapping the continuous cc mass spectrum on the physical 
charmonium resonances. 



1 Introduction 



The interplay between perturbative QCD (pQCD) and non-perturbative QCD is partic- 
ularly important in the production of charmonium states in hadronic interactions. The 
charm quark mass m,c provides a scale which is large enough to make as{m1) sufficiently 
small for the application of pQCD to describe the parton level production of a cc pair. 
The subsequent formation of a bound state charmonium state, on the other hand, is a 
non-perturbative phenomenon where important new effects enter. In the pQCD-based 
Colour Singlet Model (GSM), the predicted cross section turned out to be more than 
one order of magnitude lower than what has been observed for high-p^ charmonium at 
the Tevatron p|, |^. However, these data can be well described by the Colour Evapora- 
tion Model (GEM) |[, the Soft Colour Interaction model (SGI) § and the Colour Octet 
Model (COM) p, where the perturbatively produced cc pair in a colour octet state can 
be transformed into a colour singlet state through soft interactions which can be viewed 



as non-perturbative gluon exchange. These models make different exphcit attempts to 
describe the soft QCD djTiamics which increases the cross section for charmonium pro- 
duction by more than one order of magnitude. 

In this paper we present a detailed study of charmonium production within these 
models. For the perturbative cc production we employ exact matrix elements which 
are available up to next-to-leading order 0, ^. As an alternative we use leading order 
matrix elements complemented with parton showers giving an approximate treatment of 
arbitrarily high orders. This latter alternative is used in Monte Carlo programs which 
also include non-perturbative hadronisation resulting in the simulation of complete events. 
These event generators are phenomenologically very useful, since any physical observable 
can be extracted and compared in detail with experimental results. Our finding that this 
approach can effectively give the same result as the next-to-leading order calculation is 
therefore comforting. 

Based on these perturbative cc production mechanisms, we use the CEM and SCI 
models to form charmonium states. Here, one has to consider how the continuous cc mass 
spectrum should be mapped on the discrete spectrum of physical charmonium resonances. 
Previously used schemes have employed simple parameters that are fitted to data or factors 
from spin statistics. The latter approach is here developed further in order to introduce 
a correlation between the invariant mass of the cc pair and the masses of the different 
charmonium states. This improves the description of the observed relative rates of, e.g., 
J/ip and ip' . 

The models are tested in detail against data on J/il) production in hadronic inter- 
actions, which are available from fixed target and collider experiments covering a large 
energy range. Total cross sections and differential distributions in xp and p± can, in 
general, be well described. This gives, in particular, additional insights about important 
soft effects and should thereby help to obtain a better understanding of non-perturbative 
QCD. 

2 Theoretical description of charmonium production 

The theoretical description separates the hard and soft parts of the process based on the 
factorisation theorem in QCD. Thus, we will first consider the perturbative production at 
the parton level, then the non-perturbative formation of bound charmonium states and 
finally their combination into complete models. 

2.1 Production of cc in perturbative QCD 

Perturbative QCD should be applicable for cc production, since the charm quark mass 
TTic is large enough to make as(m^) a small expansion parameter. The leading order 
(LO) processes are gg — > cc (Fig. |l|a) and qq cc. Heavy quark production is, however, 
known to have large contributions from next-to-leading order (NLO) diagrams [0]. Virtual 
corrections to the leading order processes together with soft and coUinear gluon emissions 
give an increase of the cross section, which can be approximately described as an overall 
i^-factor multiplying the leading order cross- sect ion. In addition, some NLO tree diagrams 
with a third hard parton contribute significantly, with the most important one shown in 
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Figure 1: Illustration of cc production processes in (a) leading order and (b) next-to- 
leading order perturbative QCD. 



Fig. lb. Although this is an 0{a'l) process with a gluon splitting g —>■ cc added to the 
basic 2 2 process gg —>■ gg, it is numerically large since gg —>■ gg has a much larger 
0{al) cross section than the LO cc production processes. 

In order to treat all these contributions properly, one must use the available NLO 
matrix element calculations with explicit non-zero quark masses 0, ^. This is partic- 
ularly important at smaller energies where the correct threshold behaviour of the cross 
section is vital. We therefore use the available computer program ||^ which calculates 
total and differential cross sections based on the complete set of diagrams that contribute 
up to NLO, i.e. 0{a^), in the MS scheme. It is here possible to separate the LO 
contribution, such that the effect of the NLO corrections can be seen. The free param- 
eters of the NLO calculation are the charm quark mass rric, and the factorisation and 
renormalization scales. The latter are taken as their default options in the program, i.e. 
the transverse mass squared /x^ = + (p^) where {p\) = {plj_ +pi_L)/2 is the average 
transverse momentum squared of the charm and anti-charm quarks. The charm quark 
mass is the main parameter that affects the normalization of the cross section. The value 
rric = 1.5 GeV is found to work well in our comparisons with data below. 

The LO processes alone or with NLO virtual corrections together with soft and 
collinear emissions, cannot produce a J/ip at high-p^ • The reason is that the essen- 
tially zero p± of the initial partons cannot be conserved with the cc pair at high p± 
unless there is another hard parton to balance its p± . Of importance for high-p^ J/ 
production are therefore NLO tree diagrams having a third hard parton. The dominating 
NLO contribution is the diagram in Fig. lb. At the high energies available at colliders, 
still higher orders are expected to be important, since many gluons can be emitted and 
their virtuality need not be very large in order to allow a split into a cc pair. Higher orders 
than NLO are not available in terms of exact matrix elements, but can be approximately 
described using the parton shower (PS) approach available in some Monte Carlo event 
generators. 

As an alternative method to describe the pQCD production of cc pairs we therefore 
use the Pythia Monte Carlo, where all QCD 2 —>■ 2 processes are included with their 
corresponding matrix elements and the incoming and outgoing partons may branch as 
described by the DGLAP equations [0. A cc pair can then be produced as described by 
the LO matrix elements for qq cc and gg — cc (with explicit rric dependence) or in a 
gluon splitting g ^ cc in the parton shower. Although the latter is an approximation, it 



has the advantage of covering arbitrarily high orders and not be hmited to NLO. We use 
the defauh in Pythia for the factorisation and renormahzation scales, i.e. the transverse 
mass squared of the emerging partons in the 2 — > 2 processes (giving the same scales as in 
the NLO program). The charm mass, which affects the normalization of the cross section, 
is chosen as = 1.35 GeV to get the overall agreement with data to be shown below. 

The Monte Carlo approach gives important phenomenological advantages. Models 
for soft QCD processes can be added resulting in the generation of complete events that 
can be analysed in the same way as real events. Thus, any experimental observable can 
be extracted, including all sorts of cuts or requirements, such that a detailed comparison 
with data can be made. The NLO matrix elements, which contain corrections that cannot 
easily be interpreted in probabilistic terms, have not been implemented in a Monte Carlo 
program. Since both these approaches to describe pQCD production of cc pairs have their 
advantages and disadvantages, one should investigate to what extent they agree and can 
be used for various purposes. 

2.2 Charmonium formation in non-perturbative QCD 

The formation of bound hadron states occurs through processes with small momentum 
transfers such that the corresponding ag is large and prevents the use of perturbation 
theory. The lack of an appropriate method to calculate non-perturbative processes, forces 
us to use phenomenological models to describe the formation of charmonium states from 
perturbatively produced cc pairs. As mentioned in the Introduction, there are a few 
models that may do this based on the general idea that the colour charge of a cc pair can 
be changed from octet to singlet through soft QCD processes. A colour singlet cc pair 
with an invariant mass below the threshold for open charm {rricc < 2mj:)) will then form a 
charmonium state. In this way, the more abundantly produced colour octet cc pairs will 
contribute and increase the charmonium production by large factors. Two models of this 
kind have been investigated and are now described. 

2.2.1 The colour evaporation model 

The basic hypothesis in the CEM model ([Q and references therein) is that the colour 
exchange in the soft interactions randomise the colour charges such that no information 
remains of the colour configuration given by the preceeding hard interactions. Proba- 
bilities for colour charge states can then instead be obtained from colour SU(3) algebra, 
with the relation 3 ® 3 = 1 © 8 being applicable to a cc pair composed of a triplet and 
an antitriplet. With all colour charge states having equal weight, this implies that the 
cc will have a probability 1/9 to be in a colour singlet state and 8/9 to be in a colour 
octet state. It is then assumed that all colour singlet cc pairs with invariant mass below 
the threshold for open charm will form a charmonium state. Colour singlet states above 
this threshold, as well as the cc pairs in a colour octet state (independently of invariant 
mass) will produce open charm through the hadronisation mechanism. The cross section 
for charmonium and open charm production can then be written as 
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where nice is the invariant mass of the cc pair, rric is the charm quark mass and 2m£) 
is the DD threshold. The differential parton level cross section dacc/dmcc is the usual 
convolution of the perturbative QCD cross section with the parton density functions for 
the initial hadrons. 

This total charmonium cross section is then split on the different charmonium states 

(^i Pi (^charmonium J / '^i Vcy ~X.ci ; • • • ) ■ ('^) 

where the relative rates pi are assumed to be independent of process and energy. These 
non-perturbative parameters must in practice be determined from comparison with data. 
For J/ip, this fraction is found to be in the range pjj^ =0.4-0.5 f^. 

The CEM is one of the simplest approaches to colour neutralization, where the ef- 
fect of soft interactions is implicit in the non-perturbative factors. It gives a good phe- 
nomenological description of several production processes up to the present energy domain 
||TT| , p!2| , ^, |13[. CEM also provides an alternative to Pomeron-based models for the de- 



scription of photoproduction of J/V', D and T at HERA [14 . 



2.2.2 The soft colour interaction model 

The SCI model |T5| introduces a more explicit description of the soft interactions between 



partons and hadron remnants emerging from a hard scattering process. In these soft 
processes, the small momentum transfers are not important and can be neglected, at least 
in a first approximation. Instead, it is the colour exchanges that are important, since they 
change the colour topology of the event and thereby affect the hadronic final state. The 
model introduces an explicit mechanism where colour-anticolour, corresponding to non- 
perturbative gluons, can be exchanged between partons and remnants. This can be viewed 
as the partons interacting softly with the colour medium, or colour background field, of the 
initial hadron as they propagate through it. This should be a natural part of the process 
in which 'bare' partons are 'dressed' into non-perturbative ones and the confining colour 
flux tube between them is formed. These colour exchanges lead to different topologies 
of the confining colour force fields (strings) and thereby to different hadronic final states 
after hadronisation. 



The SCI model is implemented in the Lund Monte Carlo programs Lepto [TBI for 
deep inelastic scattering and Pythia [^] for hadron-hadron collisions. The hard parton 
level interactions are given by standard perturbative matrix elements and parton showers, 
which are not altered by softer non-perturbative effects. The probability to exchange such 
a soft gluon within a pair of partons and remnants cannot be calculated and is therefore 
taken to be a constant given by a phenomenological parameter R, which is the only free 
parameter of the model. As a result of these exchanges, the colour connections, in terms 



of strings between partons, are changed before the standard Lund model [|T^ is applied 
for hadronisation. 

The SCI model gives a novel explanation of rapidity gap events in both deep inelastic 
scattering ||15| and hard processes in pp at the Tevatron |T^, always using the same 



value R = 0.5 of the SCI parameter. Applying the same Monte Carlo implementation, 
it was found that the Tevatron data on high-p^ charmonium and bottomonium are also 
well reproduced [^. The increased production rate is here given by the possibility for a 



perturbatively produced cc pair in a colour octet state to be transformed into a singlet 
state as a result of these soft colour interactions. The mapping of cc pairs, with mass 
below the threshold for open charm production, is here made based on spin statistics 
which avoids introducing further free parameters. This results in a fraction of a specific 
quarkonium state i with total angular momentum Jj given by 

= ^ (4) 

where F = (2Jj + l)/nj corresponds to a partial width. Radially excited states are here 
suppressed through their main quantum number nj. This model was found to give a 
correct description of the different heavy quarkonium states observed at the Tevatron |^ . 



2.3 Forming the complete models 

In principle one could add the CEM or SCI models for the soft processes to any of the 
descriptions for the hard pQCD processes, i.e. the NLO matrix element program or the 
Pythia Monte Carlo generation of LO matrix elements plus parton showers. The SCI 
model is, however, constructed as an add-on to a Monte Carlo and cannot be combined 
with the NLO program. We are therefore left with three combinations to consider. 

The first model we label CEM-NLO and is the combination of the CEM model with 
the NLO program. This is done by implementing the formulas of section 2.2.1 in the 
NLO program described in section 2.1. In other words, we separate the colour singlet and 
colour octet contributions according to the CEM model and restrict the invariant mass 
region in order to produce charmonium or open charm. This gives CEM results accurate 
up to NLO in the pQCD part of the model. 

The second model is CEM-PYTHIA, where CEM has been implemented in Pythia 
version 5.7 0. The cc pairs from the perturbative phase of Pythia are then treated 
with the CEM model with respect to their colour state and projected onto charmonium 
resonances. Here, one can exploit the advantages of the Monte Carlo approach applied to 
the CEM model. 

The third model, SCI-PYTHIA, is to use the SCI model as implemented in Pythia 
5.7 0. Since this model was introduced for other purposes, namely rapidity gaps in hard 
ep and pp scattering, we are here testing its ability to describe different phenomena in a 
universal way. 

As discussed, in the NLO-based model we use rric =1.5 GeV and in the Pythia- 
based models we use rric = 1.35 GeV. These complete models also include a few other 
ingredients that need to be described. The pQCD parton level matrix elements must 
be folded with parton density distributions in the colliding hadrons in order to get the 
charm quark production cross section. To match the NLO program in CEM-NLO we 
use CTEQ5M |jl^ for protons and SMRS pOl for pions, which are parameterisations 



in NLO based on the MS scheme. The LO pQCD formalism in Pythia is matched 
with LO parameterisations. For CEM-PYTHIA, we have used CTEQ4L [^, whereas 



for SCI-PYTHIA we used CTEQ3L in order to have exactly the same SCI model as 
compared to data in previous investigations. As parton parameterisations in the pion we 



have used LO GRV-P ||23[ in both PYTHlA-based models. The uncertainty due to the 



choice of parton parameterisations has been investigated and found to be small, as will 
be discussed below. 



A further issue is the intrinsic transverse momentum k± of the partons in the colliding 
hadrons. Both in the NLO program and in Pythia, the components of a k±-vectoT 
can be chosen from a gaussian distribution. Its width is expected to be a few hundred 
MeV, corresponding to the Fermi motion within a hadron of size one fermi, but somewhat 
larger values are typically used phenomenologically. We use the width 0.6 GeV both in 
CEM-NLO and CEM-PYTHIA, whereas the slightly higher value of 0.8 GeV was used 
in SCI-PYTHIA. We note that this is in line with the new default value of 1 GeV in the 
latest Pythia version |2^. Variations of this gaussian width will be discussed below. 



The transverse momentum effect on the produced charmonium state due to the soft 
gluon exchanges that neutralize the colour of the cc pair has also been considered in the 
CEM-PYTHIA model. This was done by simply adding, to the charmonium momentum, 
a random p± chosen from a gaussian distribution with default width 0.6 GeV, as used 
for the intrinsic transverse momentum. In the SCI model, the momentum of the soft 
exchanges has so far been neglected. 

By comparisons of these three models we can separate different effects and gain under- 
standing. With GEM implemented in the NLO program and in Pythia, we can explicitly 
compare the perturbative contributions, namely NLO versus LO plus the parton shower 
approximation of higher orders. On the other hand, with SGI and GEM implemented in 
Pythia, we can explicitly compare these two non-perturbative models and see to what ex- 
tent they can account for observed soft effects and thereby provide a better understanding 
of non-perturbative QGD. 



3 Model results and experimental data 

Having discussed the details of the models in last section, we now turn to results of the 
models. Detailed comparisons between the models have been done as well as extensive 
comparison with data, both from fixed target experiments and the Tevatron collider. 



3.1 Fixed target energies 

Here we compare the models to data from different experiments using beams of protons, 
antiprotons and pions at different energies in the range 125 to 800 GeV. The targets are 
different nuclei, but the experimental results are always rescaled to represent the cross 
section per nucleon. Based on this we compare directly with our models which do not 
include any nuclear effects but treat hadron-nucleon interactions. The relatively good 
agreement between models and data, suggests that there are no strong nuclear effects 
that have been missed. 

As a first important overall test, we compare in Fig. |^ the energy dependence of the 
cross section for J/ip, ijj' and open charm obtained in the GEM-NLO model with available 
data. As can be seen, good agreement is obtained both for the overall normalization 
and the energy dependence. The GEM parameters were here chosen as pj/^ = 0.5 and 
p^r = 0.066. The former is consistent with earlier studies, whereas the latter is a new 
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Figure 2: Total cross section of J/ip-, ip' and open charm {D mesons) in proton-proton 
interactions at invariant mass \/s; curves obtained from the CEM-NLO model (colour 
evaporation model combined with NLO pQCD matrix elements) in comparison to a com- 



pilation of data |25, 26, 27, 28|. 



piece of information extracted here. The open charm cross section is obtained from all 
cc in a colour octet state plus cc in a singlet state with invariant mass above the DD 
threshold (c/. eq. (@)). This gives a result in agreement with data without additional K- 
factors. Of course, the value of the charm quark mass does affect the overall normalization 
of the charm cross section, but we note that with = 1.5 GeV the CEM-NLO model 
gives the correct result throughout the whole study presented in this paper. Although 
the other models are not explicitly shown in Fig. they can also reproduce these total 
cross sections as will be implicitly clear from the differential cross sections to which we 
now turn. 

More detailed information is available in differential cross sections. In particular, we 
will concentrate on distributions in Feynman-x, xp = p\\/p\\max in the hadronic cms, 
and transverse momentum, p±, of the produced J/ip. Fig. ^ shows these distributions for 
proton beams of different energies. As can be seen, the data are approximately reproduced, 
both in shape and normalization, by all three models. Looking into the details of the xp 
distributions, one can observe that all model curves, and in particular the CEM-NLO, fall 
less steeply than the data and therefore overshoot somewhat at large xp . The observed 
p± distribution is better reproduced, with only small differences between the models. 

Having CEM combined with different treatments of the pQCD production of cc, we 
can now investigate pQCD effects in more detail on the x^and p_L distributions. Fig. || 
illustrates this for the case of 800 GeV proton energy. 
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Figure 3: Distributions in xp and p\ of J/ip produced with proton beams of energies 800, 
530 and 300 GeV on fixed target. Data [g^ §g compared to CEM based on NLO 

pQCD matrix elements, and CEM and SCI based on LO matrix elements plus parton 
showers in the Pythia Monte Carlo. 



For the Xi? distribution in Fig. |a, the result based on the full NLO program and that 
based on LO+PS agree quite well, except for a 20% difference at small xp- The LO 
contribution in the NLO calculation is less than 50% of the total result, making the NLO 
corrections very important. In the LO+PS result, however, the PS contribution (with 
the cc pair coming from initial and final state parton showers) gives a small contribution 
which is unimportant for the overall cross section. Thus, in the PvTHlA-based models, 
the cc production is dominated by the LO 2 — >■ 2 processes having matrix elements with 
explicit charm quark mass. It is the use of a lower mass value, rric = 1.35 GeV, in 
Pythia which increases the cross section to become equal to the NLO one using the 
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Figure 4: Distributions in xp and of J / ip (in 800 GeV proton on proton as in Fig. ^) 
from variations of the pQCD treatment. GEM based on the NLO program with = 
1.5 GeV: NLO and LO matrix elements, NLO with no intrinsic k±. GEM based on 
Pythia with rric = 1.35 GeV: LO matrix elements plus parton showers (PS) and PS 
contribution shown separately. 



higher value = 1.5 GeV. We have cross-checked this within the NLO program, where 
the full result is essentially reproduced, both in shape and normalization, by the LO part 
if the lower mass rric = 1.35 GeV is used. This demonstrates that the NLO correction 
is essentially an overall i^'-factor from soft and virtual corrections. This i^'-factor can 
effectively be accounted for by using a lower charm quark mass in the LO matrix elements, 
which supports the phenomenological use of the Pythia Monte Garlo. Although, the PS 
contribution is small in the region studied here, it can become important at higher cms 
energies as we shall see in Section 3.2. 

For the understanding of the different pQGD contributions, the p± distributions in 
Fig. ^3 is even more relevant. In the LO process the cc pair cannot have a non-zero p± , 
since there is no other parton to balance its p± to obtain the essentially total zero p± 
of the initial partons. Including the intrinsic k± of the initial partons, described with a 
gaussian distribution, gives rise to a steeply falling p± distribution. The tail to large p± 
observed in data (Fig. 0) can only be reproduced when higher order pQGD processes are 
included. The NLO program gives a p± distribution with a much larger tail at large p± , 
but it is still substantially affected by the inclusion of the intrinsic k± at the limited values 
of p± accessible at fixed target energies. The p± distribution resulting from the LO+PS in 
the Pythia approach, is at high-p^ dominated by cc from gluon splittings in the partons 
showers, whereas the bulk of the cross section comes from the \o'w-p± region where the 
LO diagrams dominate. The total LO+PS result, which also includes a gaussian intrinsic 
k_i_ with the same width 0.6 GeV, agrees quite well with the NLO result. This provides 
another comforting cross-check between the NLO and Monte Garlo approaches. 

Although Fig. ^ was for the case of 800 GeV proton beam energy, the main effects 
seen and the conclusions drawn are also valid for other energies and beam particles. We 
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therefore do not show such detailed model decompositions for the following cases, but 
only compare the complete models with the data. 




Figure 5: Distributions in xp and p]_ of J/ifj produced with antiproton beams of energies 
300 and 125 GeV on fixed target. Data /|||, ^ compared to OEM based on NLO pQCD 
matrix elements, and CEM and SCI based on LO matrix elements plus parton showers in 
the Pythia Monte Carlo. 



In Fig. ^ we show the xp and p± distributions of J/ip produced with antiproton beams 
of different energies. As can be seen, all three models reproduce the general features of the 
data, both in shape and normalization, although small deviations occur. In particular, 
SCI-PYTHIA tends to overshoot the 125 GeV data. At this energy the observed p± 
distribution tends to become flatter at high-p^ , which is well described by the CEM-NLO 
model. 

Fig. 1^ shows the corresponding results for pion beams of different energies. Again, 
the models give in general a reasonable description of the data. At 300 GeV the Monte 
Carlo models describe data better, whereas at 125 GeV the CEM-NLO works somewhat 
better. 

This comparison with fixed target data on J/ip production has shown that all three 
models reproduce the main features of the data. Looking into the details of the different 
distributions, one can certainly see discrepancies such that none of the models gives an 
excellent description of all data. This is expected in view of the simplicity of the non- 
perturbative parts of the models. The observed discrepancies may therefore indicate the 
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Figure 6: Distributions in xp and p\ of J/ip produced with pion beams of energies of 
515, 300 and 125 GeV on fixed target. Data ^ |g compared to CEM based on NLO 
pQCD matrix elements, and CEM and SCI results based on LO matrix elements plus 
parton showers in the Pythia Monte Carlo. 



need for improving the CEM or SCI models. One should also remember that nuclear effects 
were neglected, which may not be appropriate for a detailed understanding of the data. 
If one lets these three models illustrate the theoretical uncertainty, such that the three 
model curves represent an uncertainty band, then there are no really significant deviations 
from the data points in Figs. H, || and ^. This shows that the essential physical effects 
in J/ip production at fixet target energies are accounted for by these model combining 
pQCD with soft colour exchanges. 

Data on ip' production provide an additional testing ground for the models, which 
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Figure 7: Distributions in xp and of ip' produced with proton and pion beams of 
energies 800 and 515 GeV on fixed target. Data |2^, ^] compared to CEM based on NLO 
pQCD matrix elements, and CEM and SCI based on LO matrix elements plus parton 
showers in the Pythia Monte Carlo. The normalization of CEM is given by p^> = 0.066 
and SCI is lowered by a factor four relative to spin statistics. 



produce all charmonium states with the same dynamics. Although there are not as much 
data on ip' as on J/ip production , one can still test the major features of the models. In 
Fig. 1^ we compare the models to data on the xp and p± distributions for ip' obtained with 
proton and pion beams. All models account quite well for the shape of the distributions. 
The proper normalization of CEM is obtained by chosing p^i = 0.066, i.e. a factor 7.6 
lower than pj/^. The spin statistics used in SCI predicts only a factor two suppression 
of ip', and must be lowered by an additional factor four in order to reproduce the data. 
Thus, although the spin statistics prescription was found to work at Tevatron energies 
1^, it fails for fixed target energies. This has prompted us to develop a more elaborate 
model for mapping cc onto different charmonium resonances which will be presented in 
section 4. 



3.2 Tevatron pp collider energy 

In order to investigate the energy dependence of the models further, we now compare 
them with data from pp collisions at ^/s = 1800 GeV in the Fermilab Tevatron. The 
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Figure 8: Distribution in transverse momentum of prompt J /if) and ip' in pp interactions 
at = 1800 GeV. Tevatron data from CDF [|| and DO [| compared to results of CEM 
added to NLO matrix elements and to the Pythia Monte Carlo. 



CDF and DO collaborations have data on the heavy quarkonium states J/ip., ip' and 
T. These can be quite well reproduced by the SCI model Here we make a similar 
comparison of the data with the results of the CEM. 

Fig. ^ shows that the CEM, both in its NLO and Monte Carlo version, can reproduce 
the data. We have here used the same p-parameters as for the fixed target case above, 
although one could optimise the overall normalization to the Tevatron data by lowering 
them by 10-20%. Such small adjustments are, however, not meaningful given the limited 
precision of the model. 

As discussed in Section 2, the LO processes for cc production are unimportant for 
high-px charmonium. This is the case here, where the higher order pQCD processes 
give the main contribution, i.e. order 0{al) tree diagrams in the NLO approach and the 
parton showers in the Monte Carlo approach. The large i^-factor associated with LO 
diagrams with soft and virtual corrections, may therefore be absent here. Indeed, the 
Monte Carlo approach does not need a lower charm quark mass to mimick i^-factors, 
but the same value rric = 1.5 GeV is here used in both CEM-PYTHIA and CEM-NLO. 
The observed shapes of these p± distributions are quite well represented by the models, 
although in particular CEM-NLO tends to give a too large cross section at the highest 
P± ■ 



The results of an extrapolation of the models to LHC energies are presented in p3| 



where also this prompt J/ip production is considered as a background for CP- violation 
studies. 
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4 Improved mapping of cc onto charmonium states 



As we have seen, the models considered give in general a good description of charmo- 
nium data both from fixed target and Tevatron energies. However, the overall energy- 
independent factors which describe the mapping of cc pair onto different charmonium 
states are not satisfactorily understood. In particular, spin statistics gives an acceptable 
ratio ip'/ ip at the Tevatron ^ , but fails at fixed target energies where an extra suppression 
by a factor four is required. In other words, this ratio has an energy dependence which has 
not been accounted for. In this section we make an attempt to understand this problem 
by developing a model which introduces a correlation between the invariant mass of the 
cc pair and the mass of the final charmonium state. 

The cc pair is produced in a pQCD process with a continuous distribution of its 
invariant mass m^c. This has then to be mapped onto the discrete distribution of physical 
charmonium states. The soft interactions that turn the pair into a colour singlet and 
form the state, may very well change its mass by a few hundred MeV, which is the typical 
scale of the soft interactions, but larger mass shifts should be suppressed. The different 
charmonium states [rjc, J/ip, Xc '^') are, however, separated in mass over a region of about 
1 GeV. It therefore seems likely that the probability to form a particular charmonium state 
will depend on the original value of nice and not just on an overall factor, such as spin 
statistics. Thus, there should be a relatively larger probability to form states that are 
nearby in mass, than those further away on the mass scale. For example, a cc with 
invariant mass just above the 2mc threshold should have a larger probability to form a 
J/ip than a ip', and a cc pair close to the open charm threshold should contribute more to 
ip' than to J/ip- Based on these considerations we have constructed the following model. 

The smearing of the cc mass due to soft interactions is described by the gaussian 

Gsme[mcc,m) = exp ( — , (5) 



2(T 



sme 



where the width agme should be a few hundred MeV. The different charmonium resonances 
i are represented by functions Fi{mi,m) around their physical masses rrij. For a cc pair 
with a given mass rricc, the probability to form a certain charmonium state is then 

X J Gsme{mcc,m)Fi{mi,m)dm SiGsme{mcc,mi) 

Vi{mcc) = ^ -, rr— — — ^ — -. r, (6) 



E j / Gsme {mcc, m)Fj {rrij , m) dm J2j Sj Gs 
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where the sum is over all charmonium resonances. The approximate result, which is used 
in our practical calculations, is obtained by letting Fi(mj, m) = Si6{m—mi), i.e. neglecting 
the width of the very narrow charmonium states and including the previously used spin 
statistical factor Sj = (2Jj + l)/nj. 

With this smearing of the cc invariant mass one must consider mass shifts across the 
threshold for open charm 2m£). Thus, a pair originally below this border may, through 
the soft interactions, acquire a mass above this threshold and hence produce open charm 
instead. Conversely, the mass could be lowered such that a pair originally above this 
threshold fluctuates downwards into the charmonium region. For a given original rricc 
value, the probability to end up in the region above DD threshold is given by the area of 
the smearing gaussian in that region, i.e. 

A{mcc) = — / exp ( — - — ^ J dm = -erfc ( — —\ , (7) 

V^Tiasme "'2mD \ ^(^sme J ^ \ y^-RO^me ) 
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Figure 9: The probability functions (1 — A{mcc))'Pi{mcc) for the different charmonium 
states as obtained in the model using the smearing width asms = 400 MeV. The resulting 
total probability for charmonium production and the remainder as open charm production 
are also shown. 

where erfc(x) is the complementary error function. The complementary area (l — A), mul- 
tiplied with the function Vi, will then give the probability to produce a given charmonium 
state. These probabilities for the different charmonium states are shown in Fig. ^. For a 
given original cc mass nice one can here see the fractional production of different states. 
The sum of all charmonium states is also shown together with the remainder giving open 
charm. One notices the non-zero contributions for charmonium also in the region above 
the DD threshold as well as some open charm production for rricc originally below this 
threshold. 

By folding these charmonium probability functions with the distribution in rricc ob- 
tained from pQCD, i.e. dacc/dnicc, the cross section for a given charmonium state is 
obtained as 

r\/s du - 

ai = / (imcc-r- ^(1 - A{mcc))Vi{mcc) . (8) 
J2mc drricc 

One can apply this mapping procedure to the GEM and SCI models and compare with 
data. In Fig. |TD|we show the ratio of ip' to J/ip production as a function of the cms energy 
of the hadron-hadron interaction. As can be seen, simple spin statistics cannot describe 
this ratio, whereas our elaborated model can. In particular, it gives a characteristic energy 
dependence of the kind indicated by the data. 

The fractions of J/ip produced directly, coming from decays of Xc states and from 
ip' are shown in Fig. |lT[ Again, the simple spin statistical factors cannot give a good 
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Figure 10: The ratio of ip' to J /ip production (times their branching ratios for decay into 
/i+yU^) in hadron-hadron interactions of cms energy y/s. Data compared to simple 
spin statistics and to our model with different gaussian smearing widths applied to CEM. 
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Figure 11: Fractions of J /if) produced directly, and coming from the decay of Xc and ip' 
states. Data ^ |^ compared to simple spin statistics and to our model with different 
gaussian smearing widths. 
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representation of all data, whereas our elaborated model gives an improved description of 
data without strong disagreements. 



The ratio Xci/Xc2 is 3 : 5 according to spin statistics. Since these Xc states are so close 
in mass, the mass correlation introduced here does not change the spin statistics result 
significantly. Furthermore, the experimental situation is very unclear for this ratio. The 
latest value from E771 is 0.53 ± 0.20 it 0.07 and they give a world average of 0.31 ± 0.14 
for fixed target experiments |3^ . In pp at the Tevatron, on the other hand, one has found 
an approximately equal production of these two states . 

We have checked that the results for J/ip presented in Section 3 do not change sig- 
nificantly if one uses this improved mapping procedure instead of spin statistics. The 
difference is less than 20%, which is within the uncertainty of the models. 



Considering all observables from Fig. 10 and 11, the preferred value for the (Jsme is 



400 MeV or lower. This fits well with the expectation for the non-perturbative dynamics 
that this model should describe. The essential new ingredient of this model, namely 
introducing a correlation between the invariant mass of the cc pair and the masses of 
the different charmonium states, is found to give an improved description of the ratios of 
different charmonium states as compared to using constant factors based on spin statistics. 
In particular, these ratios acquire an energy dependence which is particularly strong at 
low energies where threshold effects are more pronounced. A tendency for such an effect 
is observed in data, but further verification is desirable. 



5 Conclusions 

This comprehensive study of charmonium production in hadronic interactions has demon- 
strated the importance of both hard, perturbative and soft, non-perturbative QCD dy- 
namics. For the pQCD production of a cc pair, we have shown the importance of higher 
order contributions. Using a full NLO matrix element program, we find a large increase 
(factor two) of the total cross section. Since this comes mainly from soft and coUinear 
gluon emissions combined with virtual corrections, this can be effectively accounted for 
by an overall i^-factor. Because of this, the NLO results can be reproduced with the 
leading order matrix elements having a reduced charm quark mass to increase the cross 
section correspondingly. In this way the phenomenologically more useful Pythia Monte 
Carlo can be used for practical purposes. The tail of the cross section at high px of the 
cc pair is, however, dominated by higher order tree diagrams available as O^a^^) in the 
NLO program and in the parton showers of the Monte Carlo approach. 

The formation of charmonium states from a cc pair is described by non-perturbative 
dynamics. In order to understand the overall magnitude of charmonium production, one 
must here take into account that a fraction of the more abundantly produced cc pairs in a 
colour octet state are turned into a colour singlet state. This is the essence of the Colour 
Evaporation Model and the Soft Colour Interaction model considered in this study. The 
CEM uses simple SU(3) colour factors to give a random colour charge to a cc pair. For 
colour singlets with mass below the open charm threshold, different charmonium states 
are formed with probabilities given by phenomenological parameters (pj/^ etc). SCI is a 
Monte Carlo model which simulates soft gluon exchange between partons and remnants 
emerging from a hard scattering. The resulting colour singlet cc pairs (below open charm 
threshold) form different charmonium states based on spin statistics. 
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Adding CEM or SCI on the pQCD treatment provides complete models that are able 
to describe the observed charmonium production. The total cross section as well as 
the xp and p± distributions for J/ip produced with different beam particles and targets 
can be well described over the large energy range provided by fixed target experiments 
and the Tevatron collider. When using the full NLO matrix elements, there is no need to 
introduce an arbitrary X-factor since the correct normalization is obtained with the CEM 
using conventional parameter values. The SCI model, which is absolutely normalised by 
determining its single parameter to fit rapidity gap data, gives automatically the correct 
normalization. 

Looking into the details of the xp and pj_ distributions, some deviations of the model 
results from the data points can be observed. However, letting the different models 
represent the theoretical uncertainty, then there are no significant deviations. This shows 
that the essential physical effects in J/ ip production are accounted for in this kind of soft 
colour exchange models. 

The models were also tested with respect to the production of other charmonium 
states, although data is here much more scarce. The xp and p± distributions available 
for ip' can be reproduced by CEM when the free normalization parameter p^' is fitted. 
The spin statistics for the different charmonium states used in the SCI model, works 
reasonably for ip' data at the Tevatron, but gives a factor four too high ip' cross section 
at fixed target energies. Such an energy dependence of the relative rates of different 
charmonium states requires additional features of the models. 

We have made an attempt to solve this problem by a more elaborate model to map 
the cc pairs onto the physical charmonium resonances. The probability to form a certain 
charmonium state is then not only given by spin statistics, but also by a correlation 
between the cc invariant mass and the charmonium mass. Resonances with mass closer to 
the original cc are thereby favoured and fluctuations across the open charm threshold are 
considered. This improves the description of the relative rates of different charmonium 
states with an energy dependence as indicated by data. 

Our study has demonstrated that the main features of hadroproduction of charmonium 
can be described in these models combining pQCD and effects of soft colour exchanges. 
This shows, in particular, that these models for the soft QCD dynamics contain the 
essential effects and therefore improve our understanding of non-perturbative QCD. 
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